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ABSTRACT: 


Several  finite  differencing  schemes  are  used  to  compute  the  natural 
frequencies  and  mode  shapes  of  simply  supported  circular  cylindrical 
shells  in  an  attempt  to  determine  the  most  accurate  numerical  model. 
Sets  of  difference  equations  are  developed  from  the  governing  differential 
field  equations  and  by  minimizing  the  finite  difference  form  of  the 
Lagrangian  energy  function.  Both  finite  differences  and  trigonometric 
expansions  are  used  to  model  the  circumferential  behavior.  Staggered 
or  half-stations  are  used  in  addition  to  the  conventional  differencing 
schemes.   The  results  indicate  that  the  schemes  using  the  trigonometric 
expansions  are  generally  more  accurate  than  those  using  finite  differences 
for  the  circumferential  derivatives.   Furthermore,  the  conventional 
differencing  scheme  is  shown  to  be  as  accurate  as  the  half-station 
scheme  when  the  field  equation  approach  is  used  in  conjunction  with 
the  trigonometric  expansions. 
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VIBRATION  ANALYSIS  OF  CYLINDRICAL 
SHELLS  BY  SEVERAL  FINITE  DIFFERENCE  SCHEMES 

by 

Robert  E.  Ball 
Associate  Professor 
Naval  Postgraduate  School,  Monterey,  California 

INTRODUCTION 

There  are  several  large-scale  digital  computer  programs  presently  in 
operation  that  use  the  method  of  finite  differences  to  predict  the  natural 
frequencies  and  mode  shapes  of  general  shell  of  revolution  structures 
[1,2,3].  These  programs  use  a  variety  of  numerical  schemes  in  the  develop- 
ment of  the  governing  finite  difference  equations.   This  paper  uses  several 
of  these  schemes,  and  some  others,  to  compute  the  natural  frequencies  and 
modes  of  simply  supported  circular  cylindrical  shells.  Donnell's  theory 
is  employed,  and  the  unknown  dependent  variables  are  the  three  displacements, 
Two  fundamental  approaches  are  considered  in  the  development  of  the 
difference  equations;  one  is  an  energy  formulation  in  which  the  variation 
of  the  finite  difference  form  of  the  Lagrangian  energy  function  with  respect 
to  the  displacement  variables  at  each  station  is  set  equal  to  zero,  the 
other  uses  finite  difference  approximations  to  the  continuous  derivatives 
in  the  governing  field  equations.  Within  these  two  approaches,  both 
finite  differences  and  trigonometric  expansions  are  used  to  model  the 
circumferential  behavior.   In  addition,  the  three  displacements  are  either 
defined  at  the  same  stations,  as  in  the  conventional  finite  difference 
method,  or  they  are  defined  at  different,  or  staggered,  stations.  The 
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latter  technique  is  also  known  as  the  half -station  method,  and  has  "been 
used  in  the  static  analysis  of  beam-columns  [U]  and  cylindrical  shell 
roofs  [5]j  and  in  the  vibration  analysis  of  circular  rings  [6].  The 
results  presented  in  references  [U]  -  [6]  indicate  that  the  staggered 
station  scheme  yields  solutions  that  are  significantly  more  accurate  than 
the  solutions  from  the  conventional,  or  whole -station,  method.  The  goal 
of  this  paper  is  to  determine  the  most  accurate  numerical  model  from 
among  the  several  considered  for  predicting  the  natural  frequencies  and 
mode  shapes  of  simply  supported  circular  cylindrical  shells.   The 
applicability  of  that  scheme  to  the  general  shell  of  revolution  is  also 
examined . 

DIFFERENTIAL  EQUATIONS,  NATURAL  MODES, 
AND  NATURAL  FREQUENCIES 

Consider  the  thin,  elastic,  circular  cylindrical  shell  of  length  L, 
thichness  h,  and  radius  a  shown  in  Fig.  1.  The  positive  directions  of 
the  axial  coordinate  x,  the  circumferential  coordinate  8,  the  axial 
displacement  U,  the  circumferential  displacement  V,  and  the  normal  dis- 
placement W  are  as  indicated  in  Fig.  1.  Donnell's  partial  differential 
equations  for  the  vibration  analysis  of  the  cylinder  can  be  obtained  by 
minimizing  the  Lagrangian  with  respect  to  U,  V,  and  W,  and  assuming  the 
motion  to  be  harmonic  in  time  with  a  natural  frequency  Q.  This  leads  to  [7] 


> 


(1) 


u,§§  +  (-g-)  u,ee  +  (-2-)  v,§e  +  mrig  =  -  (u  u 

(— )  u,§Q  +  (-£-)  v,gg  +  v,Qe  +  w,e  =  -  u>  V 

q2  2 

-uu,g  -  v,e  -  w  -  ^  (w,§§§§  +  2w>§5ee  +  w>eeee)  =  -  *  w 

where  a  subscript  comma  denotes  differentiation,  u  =  U/a, 

v  =  v/a,  w  =  W/a,   ?  =  x/a,   or  =  h/a,  cu  =  Qa/  v  E/p(l-v  ) ,  E  is  Young's 

k 


Figure  1.  Geometry,  Coordinates  and  Displacements  of  the 
Circular  Cylindrical  Shell 


modulus,  p  is  the  mass  density,  and  v  is  Pois son's  ratio.  The  "boundary- 
conditions  for  the  simply  supported  cylinder  are 


u   =  v  =  w 


W'?§  =  ° 


§  -  0,  1 

where  I   =  L/a. 

The  coordinate  8  can  be  removed  from  equations  (l)  by  assuming 
u  =  u(|)cos  n9,  v  =  v(§)  sin  n8,  and  w  =  w(§)  cos  n0.  This  gives 


(2) 


2   /1-Un  (1+V\  2 

u,-^  -  n     (-£-)  u  +  n  K-g-)  v,     +  i?w,„  =  -  cu  u 

/1+Uv  /l-u\  2  2 

-   n  (~2~)   u'f  +  (~2~)   v>^?  "   nv~   nw=-   wv 


^  (3) 


12  V"'§S§§  "•  —  "»§g 
The  natural  modes  of  the  cylinder  that  satisfy  the  simply  supported 


-   uu,-  ••   nv  ■  ■  w  -  -75  (w,,---  -   2n2w,frfr  +  n  w)  =  -  uo  w 


boundary  conditions  (2)  are 


u  =  u  cos  (mrr§/j£)  cos  n0 

v  =  v  sin  (mrr?/i)  sin  n0  >  (h) 

w  =  w  sin  (mrr§/^)  cos  n8 
where  u,  v,  and  w  are  arbitrary  constants.  Substituting  equations  (k)   into 
equations  (l)  leads  to 
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sym. 


"12 


22 


l13 


l23 
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33 


u 


w 


=  -CJU 


u 


V 


w 


(5) 


The  coefficients  a   are  given  in  Table  I. 

pq. 

NUMERICAL  MODELS 
Consider  the  finite  difference  network  shown  in  Fig.  2a,  where  d  is 
the  nondimensional  meridional  increment  and  A  is  the  circumferential  angular 
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(b)  Model  A 


(c)  Model  B 


(d)  Model  C 


(e)  Model  D 


(f)  Model  E 


Figure  2.  Numerical  Models 


increment.  A  typical  station  is  denoted  by  i,  j,  where  i"  and  j  are  numbered 
such  that  %.   ~   id  and  9-  =  0^.  ^e   conventional  differencing  scheme  where 

J 

u,  v,  and  w  are  each  defined  at  the  same  station  is  shown  in  Fig.  2b.  A 
staggered  station  scheme  where  u,  v,  and  w  are  each  defined  at  different 
locations  on  the  shell  [5]  is  shown  in  Fig.  2c.  The  conventional  scheme 
for  the  ordinary  differential  equations  (3)  is  shown  in  Fig.  2d,  and  the 
half-station  scheme  corresponding  to  Fig.  2c  is  shown  in  Fig.  2e.  The 
staggered  scheme  used  in  reference  [2]  with  the  energy  approach  is  shown 
in  Fig.  2f.  These  five  schemes  will  be  referred  to  henceforth  as  models 
A,  B,  C,  D,  and  E  respectively. 

The  difference  approximations  to  the  partial  derivatives  for  model  A 
are  of  the  form 

(I).   .  =  (f1+i,o  -  fi-i,j>/<2d>  <6a) 

1>0 

(ft).  .  ■  <*m,i  "  2fi>d  +  fi-i,X  (6b) 

d§     1,0 


(life).   .  =  (fi+l,o+l  "  fi-l,j+l  "  fi+i,j-i  +  fi-i,;!-i)/(2a)2  (6c) 


.2 
r  ' 

60/ 
1,3 


(-JLJ: ]        =ff  +  f  +  f  +f  -2(f  +  f 

V   2    2;  L  1+1, j+i       1-1,0+1       i+1,0-1       i-i*d-i    K  i+i,J       i,d+i 

+  f .   ,     .   +  f.    .  ,)  +  kw.    .]  /d  (6e) 

1-1,0  1,0-1  IjO 


The  derivative  approximations  for  model  B  are  those  given  by  equations 
(6b),  (6d),  (6e)  and 
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(If)   =  (fi^  1  -  fi  i  i}/d  (6f) 

1  J  J 

&U  1  =  (*«.rt  " '*.*«*  "  *^.J-i  +  fi.s-i)/da  (6g) 

Corresponding  approximations  are  used  for  the  ordinary  derivatives  for 
models  C  and  D. 

Field  Equations  Approach 

The  finite  difference  form  of  equations  (l)  corresponding  to  Model  A  is 


2c 


=  -u>  u       i  =  0,1,2, ...I  ;   j  =  1,2, ...J 
-1- ?J 

(7a) 

(¥)  <Vi,j+i  -  Vi,j+i  -  ui+i,j-i  +  Vi,j-iV(^)  +  &?)  (vi+ljj  -  2vi>( 

+  Vi-l,jVd2  +   (v.?.+1  -   2vi}.   +  v^.^VA2  +   (w.5.+1  -   %iml)/(2*) 
=  "^1,3  i  =  0,1,2,...!    ;      j=l, 2, ...J 
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-K+1,0  +  6wi,j  "  Uwi-i,j  +  Vz,^/3''  +  2  [wi+i„j+i  +  wi-i,j+i  +  Vi.d-i 

+w        -  2(w      +  w.  .  ,  +  w.  ,  .  +  w.  .  n)  +  kw.     .  /(dA) 
+Wi-l,j-l    {    i+1,0    i,j+l    1-1,0    1,0-1      i>OJ 

+  (w    _  _  kw.    .  ,  +  6w.  .  -  Uw.  .  -  +  w.  .  «)/A  [■ 
k  i,j+2     1,3+1     i,0     1,0-1    i>0-2     J 

=  -m   w.  .    i  =  1,2... .1-1  ;  J  =  1,2, ...J 
i,0 

where  I  =  Z/d   and  J  =  2tt/A.   The  stations  i  =  1+1  are  ficticious  stations 
located  one  increment  off  the  ends  of  the  shell.   The  finite  difference 
form  of  the  boundary  conditions  (2)  for  model  A  is 

u_n  .  -  u.  .  .  =  v\  .  =  w.  .  =  w.  .  .  -  2w.  .  +  w.  _  .  =  0  i  =  0,1: 
l+l, j    i-l,J    IjO    ijO    1+1,0     1,0    1-1,0 

0=1,2,. .J  (7b) 

Equations  (7)  constitute  a  set  of  3(l-l)J  +  12J  homogeneous  algebraic 

2 
equations  that  define  an  eigenvalue  problem.  The  eigenvalues  are  -co  ,  and 

the  eigenvectors  are  the  natural  modes.  For  this  set  of  equations  assume 

the  modes  to  be 

u.  .  =  u  cos  (mni/l)  cos  (2nn  j/j)  (8a) 

i,0 

v.  .  =  v  sin  (mrri/l)  sin  (2nnj/j)  (8b) 

i,0 

w.  .  =  w  sin  (mrri/l)  cos  (2nirj/j)  (8c) 

1,0 

Note  that  these  modes  satisfy  the  boundary  conditions  (7b).  Furthermore, 
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substituting  equations  (8)  into  equations  (7a)  eventually  leads  to  equations 
(5)  with  coefficients  a   that  are  identical  for  every  i  and  j  defined  in 

pq 

equations  (7a) •  Hence,  equations  (8)  are  the  correct  modes  for  model  A 
and  would  be  the  eigenvectors  determined  by  a  digital  computer  except  for 
round-off  error.  The  coefficients  a   for  this  model  are  given  in  Table  I. 

pq 

A  somewhat  different  set  of  difference  equations  is  obtained  when  the 
half -station  scheme  of  model  B  is  used.  The  results  are  given  in  equations 
(Al)  in  the  Appendix.  The  natural  modes  for  this  model  are  also  given  by 
equations  (8),  but  with  i  replaced  by  i-t-§-  in  equation  (8a)  and  j  replaced 
by  j+^  in  equation  (8b).  As  in  the  previous  case,  substituting  equations 
(8)  into  equations  (Al)  leads  to  equations  (5)  for  every  i  and  j  defined 
in  equations  (Al).  The  coefficients  a   are  given  in  Table  I. 

The  same  procedure  can  be  followed  using  equations  (3),  where  the  9 
dependence  has  been  eliminated,  and  models  C  and  D.  The  difference 
equations  and  boundary  conditions  for  these  two  models  are  given  in  the 
Appendix  as  equations  (A2)  and  (A3)  respectively.  The  natural  modes 
for  both  of  these  schemes  are 

u.  =  u  cos  (mrri/l)  (9a) 

v.  =  v  sin  (mrri/l)  (9b) 

w.  =  w  sin  (mrri/l)  (9c) 

except  in  model  D,  where  i  is  replaced  by  i+^-  in  equation  (9a).  The  co- 
efficients a   for  both  models  are  given  in  Table  I. 

pq 

Equations  (A2)  are  analogous  to  the  set  of  difference  equations 
developed  by  the  computer  program  SALQRS  [3]«  The  two  sets  of  equations 
are  not  identical  because  Sanders '  shell  theory  is  used  in  SALQRS',  and 
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because  the  governing  equations  are  four  second  order  differential 

equations  in  terms  of  u,  v,  w,  and  m_,  the  nondimensional  meridional 

"bending  moment.  However,  it  can  be  shown  that  by  eliminating  the  discrete 

variable  m   from  the  difference  equations,  the  four  equation  formulation 

used  in  SALORS  reduces  to  the  identical  three  equation  formulation  given 

by  equations  (A2),  when  Donnell's  theory  is  used. 

Energy  Approach 

Sets  of  difference  equations  can  also  be  developed  by  equating  to 

zero  the  variation  of  the  finite  difference  form  of  the  Lagrangian  with 

respect  to  the  discrete  displacement  variables  and  the  Lagrangian 

multipliers .   The  Lagrangian  L  can  be  expressed  in  the  form 

8 

L=n-T+Z   A.,6,  (10) 

k=l   k  k 

where  II  is  the  shell  strain  energy,  T  is  the  shell  kinetic  energy,  X      is 

k 

a  Lagrangian  multiplier,  and  6,  is  a  constraint  condition.  Chuang  and 
Veletsos  [5]  have  shown  that  when  model  A  is  used  the  difference  equations 
developed  from  the  energy  approach  are  not  the  same  as  equations  (7a), 
which  were  developed  from  the  field  equation  approach.  However,  when 
model  B  is  used,  the  two  approaches  yield  identical  equations.  This 
section  will  consider  the  development  of  the  difference  equations  by  the 
energy  approach  using  the  models  C,  D,  and  E  where  the  derivatives  with 
respect  to  9  have  been  replaced  with  multiples  of  n. 

+  Vi 

The  shell  strain  energy  in  the  Donnell  theory  in  the  n   mode  can  be 


shown  to  be 


„ .  I^l.  f       +      d§  (lla) 

2(l-v  )   o 
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where 


o  2 

I1A  =  (u,    )     +  2v  (nv-w)  u,     +  (nv-w) 

^  »2   r,  ^2     _     2  U     2n 

12   L^W,F^  vn  ^'55       n     W   J 

%  =   ^    (V'§  "  m)2   +  1  (1"V)   ^    (w>5}' 


4*V» 

The  shell  kinetic  energy  in  the  n   mode  is  taken  as 


2       *> 


T  = 


rra  Eh 


2(l-v  )   o 
where 


TB  =  (u,t)2 


2,   J*   (TA  +  V  d?  (llb) 


and  t  =  t  v  E/p(l-v  ) /a    where  t  is  time. 

Consider  first  model  C.   The  integrals  in  equation  (ll)  can  "be 
approximated  by 


n  = 


2 

rra  Eh 


I 

1     (n.  + 


(nA  +  I^).  d  (12a) 


2(l-v2)   i=0 

I 

»-"^T7  ^  <TA  +  Vidi  ^ 

2(l-v  )   1=0 

where  (II.  +  IL).  and  (T.  +  T  ).  represent  the  strain  and  kinetic  energy 
expressed  at  the  i   station  respectively,  and  d.  =  d  except  at  i  =  0  and 
i  =  I  where  d.  =  d/2.  The  constraint  conditions  for  the  simply  supported 
shell  are  v  =  w  =  v  =  w  =  0.  Carrying  out  the  variation  of  the  Lagrangia 


Ik 


with  respect  to  the  variables  u   ,  v  ,  w   ,  u  ,  v  ,....,  v  ,  ,  w   , 
K>    ^?>  ^)  and  \,  and  assuming  harmonic  motion  leads  to  equations  (AU) 
in  the  Appendix.  A  comparison  of  equations  (AU)  with  equations  (A2) 
reveals  several  differences.  In  the  recursion  equations,  the  differences 

appear  in  the  approximation  to  the  second  derivative,  and  the  energy 

2  2  2 
approach  introduces  the  additional  term  (a   vn  d  /2U)(w.  _  -  ^w    +  6w. 

-  kw .        +  w.  p)/d  in  equation  (AUc).   Perhaps  a  more  significant  difference 
is  the  fact  that  the  five  special  equations  (Ake)   -  (AUi)  are  dissimilar  to 
the  equations  obtained  by  applying  the  boundary  conditions  (AUd)  to  equations 
(AUa)  -  (AUb)  at  i=0  and  i=l.  The  consequence  of  this  is  that  the  modes 
given  by  equations  (9)  are  not  the  eigenvectors  of  equations  (AU).  The 
fourth  boundary  condition  in  equation  (AUd)  also  eliminates  equation  (9c) 
as  the  mode  for  w. 

For  model  D,  IL  and  T  are  evaluated  at  i,  as  in  the  previous  scheme, 

.A.        A 

but  IL  and  T  are  evaluated  at  i+*r.  Hence 

I  1-1 

L    n.    a..    +  L 


n  = 


2^ 
na  Eh 


2(l-v2) 


d. 
i=0   Ai  X 


i=0 


\^ 


(13a) 


T  = 


2  , 
rra  Eh 

2(l-v2) 


I 


\.  di 

1=0    l 


1-1 

I 

i=0 


1  +  2" 


(13b) 


The  constraints  are  v  =  w  =  v  =  w  =  0.  Performing  the  variation  of  the 
Lagrangian  and  assuming  harmonic  motion  yield  a  set  of  difference  equations 
identical  to  equations  (Al). 

For  model  E? equations  (12)  are  used  in  conjunction  with  equations  (6f) 

for  (u,_).  and  (v,,).;  equation  (6a)  for  (w,  ).;  equation  (6b)  for 

I  i  51  I  1 

(w,..).;  u.    =   (u.    1    +  u.    i)/2,  and  v.    =   (v4ji  +  v.    i)/2.     The  difference 

§5    I  1  1+2  1-2  !  1+2  1-2 


15 


equations  for  this  scheme  are  given  by  equations  (A5)  in  the  Appendix. 

Note  that,  as  is  the  situation  with  model  C,  the  three  special  equations 
(A5f)  -  (A5h)  are  dissimular  to  the  equations  obtained  by  applying  the 
boundary  conditions  (A5d)  to  equations  (A5a)  -  (A5c).  Hence,  the  modes 
given  by  equations  (9)  are  not  the  eigenvectors  of  model  E. 

COMPARISON 

The  natural  frequencies  from  the  differential  equations  and  models 
A-D  have  been  computed  for  the  four  shells  a  =    .01,  1  =  2;    a  =   .01,  4=  10; 
a  -   .1,  I   =  2;  and  a  =   .1,  I  =   10;  for  several  values  of  m  and  n.  Three 
axial  increment  sizes  d  =  i/20,  JL/hO,   and  j£/100;  and  two  circumferential 
increment  sizes  A  =  2rr/l0  and  2tt/20  were  used.  Poisson's  ratio  was  .3. 
Since  there  are  three  degrees  of  freedom  per  mode,  there  are  three 
eigenvalues  for  each  value  of  m  and  n.  The  results  for  the  lowest  of  the 
three  natural  frequencies  for  the  two  shells  a   =  .01,  i  =   10,  and  a  =   .1, 
i  =  2   are  given  in  Tables  IIA  and  IIB  respectively.  Similar  results  were 
obtained  for  the  other  two  shells.  The  natural  frequencies  from  three 
dimensional  elasticity  theory  are  also  given  in  Tables  IIA.  and  IIB  [8]. 
Rather  than  present  the  natural  frequencies  of  the  numerical  models,  the 
quantity  ioo(1-uj/cd  )  is  presented  in  Tables  IIA  and  IIB  where  co  is  the 
natural  frequency  associated  with  the  differential  equations  (l).  Thus, 
the  smaller  the  number,  the  more  accurate  the  natural  frequency. 

Two  other  measures  of  the  accuracy  of  the  numerical  models  have  been 

computed.  These  measures  are 

m      n 
max    max 

ERMN  =  100   Y     V  (I-od/u,  )   /  (m    n   ) 
Li  L,  '    e'        '        max  max 

m=l    n=0 
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m       n 

max     max 


WERMN  =  100    )       )   (1-uj  /u)  )/(mn+m)   /  (m    n   ) 
L  L  mn'  e  "         '    max  max 


m=l     n=0 


where  m    =  i/d   and  n    =  rr/A.  The  two  error  measures  are  given  in 

max    '       max    ' 

Tables  IIIA  and  IIIB  for  the  three  natural  frequencies  from  each  of  the 

numerical  models,  with  n    =  5  and  10  and  m    =  20,  Uo,  and  100.  Similar 

max  max        ' 

results  were  obtained  for  the  other  two  shells. 
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CONCLUSIONS 

An  examination  of  the  results  given  in  Tables  IIA,  IIB,  IIIA,  and 
IIIB  reveals  the  following: 

(1)  Model  A  yields  extremely  inaccurate  natural  frequencies  for 
all  m  and  n  >  0  for  the  longer,  thinner  cylinder.  The  results 
are  particular  poor  for  m  -   1.  The  frequencies  are  more  accurate 
for  the  shorter,  thicker  cylinder. 

(2)  Model  B  yields  considerably  more  accurate  natural  frequencies 
than  model  A  for  the  longer,  thinner  shell,  and  about  equally 
accurate  frequencies  for  the  shorter  thicker  cylinder,  except 

for  m  =  1  where  model  B  is  the  more  accurate.  In  general,  model 
B  is  the  more  accurate  model  for  all  four  shells. 

(3)  Models  C  and  D  yield  essentially  the  same  natural  frequencies 
for  all  modes  for  all  four  shells .  In  some  instances ,  model  C 

is  slightly  more  accurate,  and  in  others,  model  D  is  the  more 
accurate.  In  general,  models  C  and  D  yield  more  accurate 
frequencies  than  models  A  and  B. 

Hence,  model  B  is  preferred  over  model  A,  and  models  C  and  D  are 
preferred  over  models  A  and  B,  provided  the  field  equation  approach 
is  used  to  develop  the  difference  equations  for  model  C.  The  conclusion 
that  model  B  is  more  accurate  than  model  A  is  in  agreement  with  the 
results  presented  in  reference  (5).  The  conclusion  that  models  C 
and  D  are  of  equal  accuracy  is  somewhat  surprising,  in  view  of  the  fact 
that  model  D  yields  identical  difference  equations  from  both  the  field 
equation  and  the  energy  approach. 
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The  staggered  station  schemes  of  models  B  and  D  appear  to  be 
ideally  suited  to  Donnell's  theory  with  the  boundary  conditions  given 
by  equation  (2)  because  only  even  or  odd  ordered  derivatives  of  the 
dependent  variables  appear  in  any  one  field  or  boundary  equation. 
General  shells  of  revolution  with  meridional  curvature  do  not  have  the 
consistant  even  and  odd  ordering  of  derivatives.  Hence,  the  staggered 
station  schemes  do  not  appear  to  be  more  suitable  for  the  general 
shell  of  revolution  than  for  the  cylinder.  On  the  other  hand,  the 
meridional  curvature  may  have  some  influence. 
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APPENDIX 
Difference  Equations  for  Model  B 

The  equations  at  a  typical  station  are 

p  Id 

+  ui4,d-i)/A2  +  (^>  (vi+i,d4  -  vi,i4  •  vi+i,a-i  +  vu-i)/(aA) 

+  u(w-i+i    -i   "  wi    nVd  =  -   ">  u.    i    j  i  =  0,1,2.. .1-1   ;     j  =  1,2,  ...J        (Aia) 


2  2 

-   2v.     .    i   +  v.    n     .   i)/d     +   (v.     .  3    .   2v.     .,i   +  v.     .   i)/A 

2 
+  (w.    .   .    -  w.     .)/A  =  -  (i)  v.    .  i  i  =  1,2,... 1-1  ;       j   =  1,2,... J       (aj+j) 

1,<J+-L  -LjJ  •J-,J"'"2 


2 

-  *(ui4.3  -  Vi.jVa  -  (vi;j4  -  v1)3.#)/A  -  .±>J  -  2  {(wi+2!j 

-   W.    .,     .   +  6w  .     .   -W.    n    .   +  w.    0    .)/d     +2       w.    ,     .   ,    +  w.    _    .    ,    +  w.    .     .    .. 
i+l,j  i,J  i-1jJ  i-2, Jy/  L    l+l, j+1         i-l,J+l         i+l, J-l 


+  w 


Wi)J+2  "  Uwi,d+1  +  6wi;j  "  N.J-1  +  "i,i.ZV*k} 

=  -  (l>  w.;j     i  =  1,2,.. .1-1  ;    j   =1,2,  ...J  (Ale) 
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The  boundary  conditions  are 

u.    i    .   -   u.    i    .   =  v.     .   i   =  w.     .   =  w.    .    .   -2w.    .   +  w.    .,    .   =  0 
1+2JJ  1-2>0  ijJ+2  i»J  1+1jJ  i»J  i-l, J 


i  =  0,   I   ;   j  =  1,2, ...J  (Aid) 

Difference  Equations   for  Model  C 

The  equations  at  a  typical  station  are 

(ui+i  -2ui +  viV*2  - n2  <¥>  ui +  n  <¥>  (vi+i  -  vi.1)/(^) 

+  u(w1+1  -  w±_1)/(2d)     =     -A.  1-0,1,2,..-!     ;n  =  0,1,2,...      (A2a) 

-n  (^)   (u.+1  -   u._1)/(2d)  +  (2$)   (v.+1  -2v.    +  v.^Vd2     -  n2v. 

"   nwi  =  "   *\  i  =  0,1, 2,... I      ;   n  =  0,1,2,...     (A2b) 

2 
"   U(ui+1  -  VlV(2d)  -nv.   -  w.    -  £-       (w.+2  .   Wi+i  +  6Wi  .  ^  +  V2) 

/d       -2n     (w         -  2w.    +  w.    .  )/d2  +  nV 
i+i  i         l-l  '  i 


2 

=  -co  w. 

1 


i  =  1,2,.. .1-1      ;      n  =  0,1,2,...   (A2c) 


The  boundary  conditions  are 

Vi  "  ui-i  =  vi  =  wi  =  Vi  -  2wi  +  Vi  -  ° 


i  =  0,1  (A2d) 


Difference  Equations   for  Model  D 

The  equations   at   a  typical  station  are 


(\4  -  ^4  +  V^/d  -  n   <¥>  ui4  +  n  (^}  (vi+i  -vi)/d 

(A3a) 
+  v   (»         -  WjL)/d       =  -   uj  ui+i  i  =  0,1,2,...I-1   ;   n  =  0,1,2, 


-n  (*£)   (u.4  -   ui4)/d  +  {*£)    (v.+1  -   2v.    +  v.^Vd2       -n2  v. 

(A3b) 

2 
-   nw±  =  -  u)  vi  i  =  l,2,...l-l  ;   n  =  0,1,2,.. 


2     _ 
-   u  (ui4  "  ui4}/d  -   nv.    -  w.    -  -g     [(w1+2  -  %ri+1  *  6w. 


"   H-l  +  wi-2^d     "  2n2  (wi+l  "  ^i  +  wi-lVd2  +  n\] 


(A3c) 


2 
=  -  w  Wi  i  =  1,2,.. .1-1  ;   n  =  0,1,2. 


The  boundary  conditions  are 

u.  i  -  u.  i  =  v.  =  w.  =  w.  ,  -2w.  +  w.  ..  =  0 
1+2    1_2    1    x    1+1    x    1-1 


i  =  0,1  (A3d) 
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Difference  equations  for  Model  C,  Energy  Approach 

The  equations  at  a  typical  station  are 
(u.+2  -  2u.  +  u._2)/(2d)2  -  n2  (2gH)  u.  +  n  (^)  (v.+1  -  r±   ,)/(2d) 


+  u(wi+1  -  wi_1)/(2d)  =  -  oj2Ui     i  =  2, 3,... 1-2  ;  n  =  0,1,2,...  (A^a) 


"  n  (1T]    (Vl  "  VlV(2*)  +  (i6S)  (Yi+2  "  2vi  +  vi.2V(2d)2  -  n2v. 


nwi     =  -  uu  Vi  i  =  2, 3,... 1-2  ;  n  =  0,1,2,...  (Akb) 


-  v   (u.  n  -  u.  ,)/(2d)  -nv.  -w.  -  T5   (l  +  un  d  /2)(w.  0  -Uw.  .  +  6w.  - 
v  l+l    1-1"  v       l   l   12  Lv         '  /v  i+2    l+l     i 

Wi_1  +  wi_2)/d4-  2n  (wi+2  -2wi  +  wi_2)/(2d)  +  n  v±   |   =   -  cu  Wi 


=  2,3, ...1-2  ;   n  =  0,1,2,...  (A4c) 
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The  equations  corresponding  to  the  boundary  conditions 

at  the  initial  end  are 

2  2 
u,  -  u,  =  vQ  =  Wq  =  w,  +  w  --  (l-u)n  d  (w,  -  w_1)/2  =  0  (A^d) 


There  are  five  special  equations 

(u2  -  uQ)/(2d)2  +  unV;L/(2d)  +  uw1/(2d)   =  -  umi^       (A^) 


(vl  "  v_iV(2d)  "  n^0  =0  (A4f) 


(u3  -  u1)/(2d)2  -n2  (^)  ux  +  n  (^)  v2/(2d) 


+  uw2/(2d)  =  -  ou2^  (A^g) 


-  n  (^)  (u2  -  uQ)  -  n  (%£)  uQ/(2d)  +  (^)  (v3  -  v1)/(2d) 


2  2 

+  n  v  -  nw..  =  -  uj  v  (A^h) 

-u  (u2  -   uQ)/(2d)  -nv1  -w1  -  ^    {  [w3  -  k22  +  (11/2  -  Y)  wj  /dU 
+  (un2d2/2)     [w3  -   W2  +  (13/2  -   v)  wx  J/d^  -  2n2      [w3 

-    (3/2  +  Y)  wx  ]/(2d)2  +  n^  }  =  -  A^  (Ml) 

where 

Y  =  ■§    [2  -    (1-u)   n2d2]   /  [2  +  (1-u)   n2d2  ] 
Similar  equations  apply  at  the  final  end  of  the  cylinder. 
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Difference  Equations  for  Model  E,  Energy  Approach 
The  equations  at  a  typical  station  are 


(\4  -  2ui4  +  ui4)/d2  "  °2  (tt>  (ui4 +  2ui4  +  va'A 


r%  A-rg-  X-2 


■  1+v 


+  n  (±g.)   (vi+_3  -   v.  ^)/(2d)   +  u   (wi+1  -   Wi)/d  = 


-   u,2   (ui+3  +  2u.4  +  ui4)/lt 


i  =  1,2,.. .1-2  n=  0,1,2,...       (A5a) 


-n  C^>  (u1+|  -  ui4)/(2d)  +  (^)  (vi+|  -  2v.4  +  vi4)/d2 

14-o  1+3"  l-"5"  l+±  1  1-Tp  X-r*  1-p 


i  =  1,2, ...1-2  n=  0,1,2,...      (A5b) 

2 

"u(ui4  "  Ui-i}/d  "   n(vi4  +  vi-i}/2     -*i  -  ig  (l^n2d2/2)    (w.+2 

-W.+1  +  6Wj_  -   Wi_1  +  w._2)/dlf  -2n2   (wi+2  -   2w.  +  wi_2)/(2d)2 


li  p 

+  nw.  =_uuw.  i  =  2,3, ...1-2  n  =  0,1,2,...        (A5c) 
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The  equations  corresponding  to  the  boundary  conditions  at  the 
initial  end  are 

2   2 
Vi  +  v  i  =  wQ  =  w-   +  w  -  -    (l-u)   n  d     (w,    -  w  ,  )/2     =     0  (A5d) 


(ux   -   u  x)/c?  -n2   (&2)   (ui  +  u  i)/l*  +  n(i^)  vi/d 

2  "2  ^  2  "2  ^2 


=  -w2  (ui  +  u  i)/if  (A5e) 

2  "2 


There  are  three  special  equations 

(u3  -  -Jui     +  Ju  i  )/d     -n     (-^)    (u^  +   |  u i  +  §  u  i)/U 

2  2  "2  ^2^2  "2 


+n  (^)   (v3  +  V!)/(2d)   -   n  (^)  vi/(2d)  +  u  w ,/d  = 

^22  ^2  J- 


2     ,     „  3 


-  ur  (u3  +  fui   +  Ju  i)  (A5f ) 

2  2  "2 


_n(l±H)   (u3  .   u  i)/(2d)   +  un  (ui   -   u  i)/(2d)   +  (i^) 

d.  2  -2  2  "2  c 


(y3    -  yi  )/d2  -n2   (vj  +  vi  )/4  +  nw-,/2     =     -uu2(y3  +  v±)/k  (A5g) 

2  2  2  2  2  2 


-u(u3  -  ui)/d  -n(v3  +  vi)/2  -  wx  -  ~   {  [w3  -   4w2  +  (ll/2  -  y)  ^  ]    /d4 

+  (»n2d2/2)    |"  w3  -   W2  +   (13/2  -   y)   "J/cl4  -   2n2[  w3  -    (3/2  +  y)  wj  /(2d)! 

if       \  2 

+  n     w-J     =  -   u>     w. 

Similar  equations  apply  at  the  final  end  of  the  cylinder. 
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